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Abstract. We derive the numerical schemes for the strong order integration of the set of the 
stochastic differential equations (SDEs) corresponding to the non-stationary Parker transport 
equation (PTE). PTE is 5-dimensional (3 spatial coordinates, particles energy and time) Fokker- 
Planck type equation describing the non-stationary the galactic cosmic ray (GCR) particles 
transport in the heliosphere. We present the formulas for the numerical solution of the obtained 
set of SDEs driven by a Wiener process in the case of the full three-dimensional diffusion tensor. 
We introduce the solution applying the strong order Euler-Maruyama, Milstein and stochastic 
Runge-Kutta methods. We discuss the advantages and disadvantages of the presented numerical 
methods in the context of increasing the accuracy of the solution of the PTE. 


1. Introduction 

We employ the stochastic methodology to model the galactic cosmic rays (GCR) transport in 
the heliosphere. Modulation of the GCR is a result of an action of four primary processes: 
convection by the solar wind, diffusion on irregularities of HMF, particles drifts in the non- 
uniform magnetic field and adiabatic cooling. Transport of the GCR particles in the heliosphere 
can be described by the Parker transport equation [T]: 

(K^J ■ V/) - (Fd + c7) • V/ + |(V ■ U)^, (1) 

where / = f{f,R,t) is an omnidirectional distribution function of three spatial coordinates 
r = {r,9,ip), particles rigidity R and time t] U is the solar wind velocity, Vd the drift velocity, 
and Kfj is the symmetric part of the diffusion tensor of the GCR particles. 

This paper is an extension of our previous results presented in [2]. We have presented [2] that 
the GCR transport can be effectively modelled based on the solution of the set of stochastic 
differential equations (SDEs) corresponding to the PTE ([1]). Firstly, the PTE must be brought 
to the form of the backward Eokker-Planck equation (e.g. m), and then the corresponding SDEs 
must be solved numerically. In [am the Euler-Maruyama method was applied. In this paper, 
we increase the accuracy of the SDEs solution by applying the higher order methods i.e. Milstein 
and stochastic Runge-Kutta. 


2. Numerical methods to solve the stochastic differential equations 

The Euler-Maruyama method of solution of the SDEs is unconditionally stable in the higher 
dimensions and doesn’t rely on the numerical grid. All integration method involves the statistical 
error, which can be reduced by increasing the number of simulated pseudoparticles and applying 
the higher order integration method e.g. Milstein’s method or Runge-Kutta method e.g.[5]. 
These three methods improve the approximation of the solution by applying higher order 
extension of the SDE solution in the Ito - Taylor series e.g.[U]. The Euler-Maruyama method 
has order of convergence 7 = 0.5, addition just one more term from the Ito - Taylor expansion 
the order of convergence increases for the Milstein method up to 7 = 1. In turn, the stochastic 
Runge—Kutta method requires to generate a new random variable Z{t) resulting in strong order 
of convergence with 7 = 1.5. 

Let’s consider a basic SDE: 


dXt = f{Xt)dt + g{Xt)dWt, (2) 

depending on the applied method the numerical approximation of the solution can be obtained 
by following formulas : 

Aj_|_i = Xj + f ■ dt + g ■ dWt + 2 ^ ’ 9 — dt) -|-$ (3) 

Euler-Maruyama 

" -V-" 

Milstein 

" -V-" 

StochasticRunge—Kutta 

where $ = f'-g-dZj+\{f -f' + \g'^-f'')dt^ + {f-g' + \g'^-g'')-{dWj-dt-dZj) + \g{g-g"+g''^){\dWf-dt)dWj 
and dZi = ^dt{dWi + and dVi is an additional Wiener process e.g. [5] . 


3. Stochastic differential equation corresponding to the Parker transport equation 

The PTE (EqlT]) in the 3-D heliocentric spherical coordinate system {r, 9, tp) written as time-backward 
FPE diffusion equation has the form: 


, d^f . dV ^ dV ^ d^f , d'^f , d^f 

dt ~ dr 2 dv ?2 + ^4 + As + As 


A 9f df df A 


the coefficients Ai,...,Aio are presented in detail in [2]. Depending on the choice of the numerical 
approximation the corresponding to EqjJ] set of SDEs has a form: 
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dR = Aio • dt. 
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Figure 1. Trajectories of the pseudoparticles with rigidity 10 GV for the A>0 solar magnetic cycle 
obtained applying the Euler-Maruyama, Milstein and Stochastic Runge-Kutta method. The specific 
colors highlight the trajectories of the sample pseudoparticles, based on the same Wiener process, traced 
backward in time from the heliosphere boundary until they reach the position r = lAU, 9 = 90°, 

ip = 180°. 


The Bij, {i,j = r,9,ip), is the lower triangular matrix presented in [2], and the coefficients 
<hi, <f> 2 , $3 and $4 have the form: 

chi = Bn • dZr^ + + lB‘l^^)dt^ + {A^^ + \Bl^^){dWr • dt - dZr) + 

iBniBu^ + {^?)C 3 dW^ - dt)dWr-, 

Ch2 = Bn • dZr^ + + ^Bl^)dt^ + {As^ + ^Bl^){dWr • dt - dZ^) + 

\B2i{B2i^ + {^?){\dW^-dt)dWr + B22-dZe^ + }^{A^^ + }^Bl^^)de + {As^ + 
\Bl 2 ^){dWe • dt - dZg) + \B22{B22^ + {^?){\dWl - dt)dWg; 

^3 = - dt) + lB,2^{dWi - dt) + iBss^idW^ - dt); 


Ch4 = Bn ■ dZr^ + l(Ag^ + ^^Bl^)dt^ + (Agff + \Bl,^){dWr • dt - dZr) + 
^^Bn{Bn^ + {^)^){yW?-dt)dWr + Bs2-dZg^ + ^,{A,^ + ^,Bl^)dt^ + {A,^ + 
hBi,^){dWg-dt-dZg) + ^^Bn{Bs2^H^?)C3dW^-dt)dWe+Bs3-dZ^^ + ^^{A,^ + 
lBl^)dt^+iA,^+^Bl^)idW^-dt-dZ^)+^^BM,^+i^)^)ildW^-dt)dW^. 


4. Results and Summary 

We performed the simulation applying all three methods given by Eqsi5l The trajectory of 
n = 3000 pseudoparticles was traced backward in time in the spherical heliocentric coordinate 
system. The pseudoparticles were initialized in the point representing the Earth’s orbit (i.e. 
r = lAU, 6 = 90°, ip = 180° ) and traced backward in time until crossing the heliosphere 
boundary assumed at 100 AU (see Eig. 1 in [2]). The value of the particle distribution function 
f{f,R) for the starting point was be found as an average of fiisiR) value for pseudoparticles 























































Figure 2. (a)Simulated galactic protons rigidity spectra and the histograms of (b) the particles rigidity 
and (c) the particles exit time for the pseudoparticles initialized R=10 GV from position r = lAU, 
0 = 90°, (p = 180° for the A>0 solar magnetic cycle obtained applying the Euler-Maruyama, Milstein 
and Stochastic Runge-Kutta method. 


characteristics at the entry positions, f{f,R) = jfJ2n=ifLis{R), where fiisiR) is the cosmic 
ray local interstellar spectrum taken as in [8] for rigidity R of the particle at the entrance 
point. 

In the case when we know the analytical solution it is easy to compare the efficiency of each 
numerical method. However, here we do not have such a possibility. Thus, we have compared 
the simulated galactic protons rigidity spectra for A>0 solar magnetic cycle. Fig. [2^ presents 
that all applied numerical schemes resulted in the same spectra, slight differences are seen for 
the lower rigidity particles. The obtained spectrum is in agrement with results presented in [7j- 
Fig. [1] presents the trajectories of the pseudoparticles with rigidity 10 GV for the A>0 solar 
magnetic cycle with respect to all coordinates. To allow the reliable comparison all simulations 
were based on the same Wiener process. One can see that there arises some subtle differences: 
the trajectories of pseudoparticles vs. the heliolongitude are the widest for Euler-Maruyama and 
narrow for the Runge-Kutta method. At the same time the opposite distribution is observed 
vs. the heliolatitude. Additionally, we have estimated the histograms of the particles rigidity 
and the particles exit time for all three methods. Fig. [2jrc show that the time of traveling the 
pseudoparticles throughout the heliosphere is the longest applying the stochastic Runge-Kutta, 
at the same time particles loses less energy during travel. 

The obtained results suggest that all applied methods give a reliable results. Performed 
tests proved that the Milstein and stochastic Runge-Kutta methods are more stable and return 
the same values of differential spectra when we decrease the number of simulated particles by 
factor of three. However, the slight differences between the pseudoparticles trajectories and its 
characteristics (rigidity and exit time) given by the tested numerical schemes must be analyzed 
in more detail. 
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